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We show that the hnear dynamics of cosmological perturbations can be described by coupled wave 
equations, allowing their efficient numerical and, in certain limits, analytical integration directly in 
position space. The linear evolution of any perturbation can then be analyzed with the Green's 
function method. Prior to hydrogen recombination, assuming tight coupling between photons and 
baryons, neglecting neutrino perturbations, and taking isentropic (adiabatic) initial conditions, the 
obtained Green's functions for all metric, density, and velocity perturbations vanish beyond the 
acoustic horizon. A localized primordial cosmological perturbation expands as an acoustic wave 
of photon-baryon density perturbation with narrow spikes at its acoustic wavefronts. These spikes 
provide one of the main contributions to the cosmic microwave background radiation anisotropy 
on all experimentally accessible scales. The gravitational interaction between cold dark matter 
and baryons causes a dip in the observed temperature of the radiation at the center of the initial 
perturbation. We first model the radiation by a perfect fluid and then extend our analysis to account 
for finite photon mean free path. The resulting diffusive corrections smear the sharp features in the 
photon and baryon density Green's functions over the scale of Silk damping. 



I. INTRODUCTION 

The nearly perfect black body spectrum and isotropy 
of the cosmic microwave background radiation indicates 
that the universe at large redshift was highly uniform 
and in thermal equilibrium, with fluctuations in temper- 
ature of only about 1 part in 10^ on observable scales. 
Under these conditions, the dynamics of matter, radia- 
tion, and gravity is described accurately by linearizing 
the governing equations about their spatially homoge- 
neous solutions representing an unperturbed expanding 
background. In place of the strongly nonlinear fluid and 
Einstein equations, we have a system of coupled linear 
partial differential equations. The linear approximation 
continues to apply much later on large scales even af- 
ter nonlinear structures such as stars and galaxies have 
formed on smaller scales. 

The cosmological perturbations can be described by a 
set of classical fields, e.g. the gravitational potential 0. 
After the perturbations are created in the early universe, 
each field may be expanded over a convenient set of basis 
functions 

</)(r,r,) = ^0fc/fc(r) , (1) 

k 

where the (pk are expansion coefficients at some initial 
time r,;. Under linear evolution at a later time r 

0(r,T) =^^fcFfc(r,T) , (2) 

k 



where (r, r) is the solution to the linearized equation 
for (f) that satisfies the initial conditions Ffc(r, t^) = /fc(r). 

The range of possible basis functions /fc(r) is enor- 
mous. Since the pioneering work by Lifshitz, Ref. [0, 
nearly all the cosmologicalperturbation theory calcula- 
tions, e.g. Refs. [|], ||, ^, || ||, 0|, have used harmonic 
plane waves or their generalizations in curved spaces. 
There is a good reason for this: because the dynami- 
cal equations are translationally invariant, each harmonic 
mode evolves unchanged aside from a time-dependent 
multiplicative factor Tkir) called the transfer function: 
Ffc(r,T) — Tk{T) fk{r). This separation of variables al- 
lows the partial differential equations to be reduced to a 
set of ordinary differential equations that are straightfor- 
ward to solve numerically. 

However attractive this reduction appears, it by no 
means implies that cosmological perturbation theory re- 
duces to the simplest form in Fourier space. The plane 
wave expansion may be thought of as a localized (Dirac 
delta) basis in Fourier space. From this perspective, it 
is not unreasonable to consider a localized basis in real 
space. In this case, the function Ffc(r,r) in eq. (||) no 
longer factors into a simple product. The lack of separa- 
tion of variables would seem to imply that perturbation 
theory is more difficult in position space. In fact, we will 
show that it is not only simple in the linear regime, but 
the dynamics is clearer and more intuitive in real space. 

Suppose, for example, that we consider the evolution 
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of a perturbation originating from a point disturbance 



(3) 



where S§\r - Tq) = 6b (x - Xo)6Biy - 2/o)^d(2 - zq) is 
the product of Dirac delta functions. What will such an 
initial perturbation evolve to by time r? The correspond- 
ing function Ffc(r, r) = G{r — ro,T) is called a Green's 
function. To see how the Green's function language can 
be simpler than the transfer function approach, we no- 
tice that, as shown in Ref. prior to recombination 
the primordial isentropic (adiabatic) perturbation satis- 
fying the initial conditions (^) expands in the photon- 
baryon plasma as a spherical acoustic wave with sound 
speed Cg. During the radiation era, Cg = c/VS and the 
gravitational potential is simply given by the properly 
normalized Heaviside step function: 



r, T 



47r(csT)" 



r < CsT 
r > c,T 



(4) 



Here, r is the comoving distance from the initial point 
perturbation, and r is conformal time related to the 
proper time thy dr = dt/a{t), where a{t) is the cosmolog- 
ical scale factor. Later we will also use Green's functions 
produced by initial disturbances on two-dimensional 
planes in space. 

The position space representation is formally equiva- 
lent to the Fourier space representation. It can be used 
to describe the dynamics of cosmological perturbations 
regardless of their origin and statistical properties. The 
Green's function method can be applied just as well to 
the evolution of nearly scale-invariant perturbations gen- 
erated by inflation, Refs. I^], ||], as it can be apphed 
to local perturbations produced by topological defects, 
Refs. 00. 

So what new can be learned from this approach? 
Green's functions contain the same information as trans- 
fer functions, but that information is packaged differ- 
ently. The Green's function approach reveals new sides 
to cosmological dynamics that are of both phenomeno- 
logical and theoretical interest. 

Phenomenologically, Green's functions are often char- 
acterized by localized features such as the acoustic wave- 
front. Through the uncertainty relation AxAk ^ 1, 
localization in position space results in features being 
spread over a broad range of wavenumbers in Fourier 
space. An example of this is the acoustic peaks in the 
cosmic microwave background (CMB) power spectrum 
Ci. In Ref. 1^ we have shown that the position-space 
analogue of Ci, the angular correlation function C{9), has 
localized features instead of acoustic oscillations. These 
features offer an alternative signature for experimental- 
ists to measure. 

Theoretically, acoustic and transfer processes are made 
explicit in position space. This offers new methods for 
solving the evolution equations or leads to substantially 
simpler equations and solutions in certain cases. An ex- 
ample is the spherical wave solution for the radiation era 



given by eq. (Q). The acoustic wave is manifest, and 
this suggests that the position space view may provide 
a more direct understanding of the dependence of CMB 
anisotropy patterns on the underlying cosmological pa- 
rameters. 

Equations of perturbation dynamics can in principle 
be numerically solved faster in position space than the 
equivalent equations in Fourier space for same desired ac- 
curacy. This is because the Green's functions are mono- 
tonic and limited in their spatial extent by the acoustic 
horizon for perfect fluids, while the Fourier transfer func- 
tions oscillate in both the wavenumber and time coordi- 
nates requiring a larger number of sample points for their 
accurate representation. 

This paper describes cosmological linear perturbation 
theory in position space using the Green's function ap- 
proach. The discussion is primarily focused on the pe- 
riod of cosmological evolution prior to hydrogen recom- 
bination and radiation decoupling from baryonic matter. 
In the current paper we derive evolution equations that 
are convenient for a position space analysis and consider 
their Green's function solutions. In a later paper we will 
show how these results may be used to describe CMB 
anisotropy. Sec. || presents the dynamical equations de- 
scribing coupled perturbations in the metric and in radi- 
ation and matter modeled by locally isotropic fluids. In 
Sec. Ill wc specify the initial conditions for the Green's 
functions and give explicit Green's function solutions in 
the radiation epoch. In Sec. IV we discuss general prop- 
erties of the Green's functions and their numerical inte- 
gration in the fluid model, and we analyze the results of 
numerical integration up to the time of hydrogen recom- 
bination. Sec. ^ shows how the position-space descrip- 
tion of perturbations can be implemented when the fluid 
equations of Sec. || are replaced by the full equations of 
Boltzmann phase space dynamics, which are summarized 
in the Appendix. Then in Sec. |V| we find the leading 
corrections to the fluid approximation results. A brief 
summary is given in the Conclusion. 

We end this introduction with a summary of our con- 
ventions and notations. Throughout the paper, Greek 
indices range from to 3 and label components of space- 
time vectors. Components of spatial 3- vectors carry 
Latin indices ranging from 1 to 3; if the indices are omit- 
ted, 3-vectors are typed in bold. The speed of light is 
c = 1. The 2TT factors in the Fourier transforms always 
appear in the denominator of the momentum integral, as 
in the equation 



(r) 



(2^ 



(5) 



We consider only the scalar perturbation mode, the one 
involving radiation and matter density perturbations, 
and work in the conformal Newtonian gauge for the per- 
turbed Robertson- Walker metric, Refs. |Q, |^. In this 
gauge 



ds-" 



a^ir) [-(1 + 2(t))dT^ + (1 - 2tl))dr^ 



(6) 
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Symbol 


Meaning Defi 


tiing Equation 


a 


Scale factor relative to the present 




y 


Scale factor relative to the radiation-matter density equality 




T 


Conformal time 




Te 


Characteristic r of the radiation-matter density equality 


w 


r 


Comoving 3-space coordinate 




allJJclbCliJJL 1 O 1 


Q Vl i^T"l^^Ql 1*^ 1 t ^ T"£i^^Tl 'c fllTl^^^"l/^Tl 

DpiiCliCdl j Vjiceii ILlllCLiOli 


\iS.) 


superscript (1) 


Plane-parallel (ID) Green's function 


(1) 


X 


Spatial coordinate of a plane-parallel Green's function 


(3£) 


superscript (T) 


Fourier space transfer function 


m 


k 

ft 


Comoving wave vector 


m 


0,V' 


|k| 

Metric perturbations 


- 

® 


subscript 7 


Photons 




subscript b 


Baryons 


- 


subscript r 


Radiation fluid (coupled photons and baryons) 


- 


subscript d 


Difference in photon and baryon perturbations 


- 


subscript 1/ 


Massless neutrinos 


- 


subscript c 


Cold dark matter (CDM) 


- 


subscript rel 


Background of relativistic species (7 plus u) 


- 


subscript m 


Background of non-relativistic matter (baryons plus CDM) 


- 


6 


Energy density enhancement / 


(lOa) 




If applied to a variable, the perturbation of that variable 


- 


Su{x) 


Dirac delta function 


- 


u 


Peculiar velocity potential 


(lOH) 


Vi 


Peculiar 3-velocity 


jioil) 


5p 


Pressure perturbation 


(|i6c|,|A2|) 


w 


Shear stress potential 


(|A2|) 


a 


Entropy perturbation potential 


(||® 


f 


Phase space distribution 


(Ey 


F 


Energy averaged perturbation of / 




fi 


Potentials for angular moments of F 


(|A1C|) 




Baryonic fraction of nonrelativistic matter (pt/pm) 




Cs 


Speed of sound in the photon-baryon plasma 




Cw 


1/2 

"Isentropic sound speed" {dp/dp)J^^^^ 


(} 


A and B 


1/(3c55j) and l/(3Cs) respectively 


(2£) 


Sir) 


Radius of acoustic horizon 


(1£) 


Tc 


Mean conformal time of a photon collisionless flight 


(75) 



TABLE I: Frequently used notations. 



where dr^ is the three-metric of a Robertson- Walker 
space. For perfect fluids, there is a single gravitational 
potential 4> — ip, but in general two distinct potentials 
are required. (Note that the perturbation (j) of the met- 
ric (H) is called i/', and V' is 0, in Refs. |§, Our present 
choice agrees with Ref. [||.) Other frequently used vari- 
ables and notations are summarized in Table 1. They 
will be introduced systematically in what follows. 



II. COSMOLOGICAL DYNAMICS IN THE 
FLUID APPROXIMATION 

A. The model 

In this and the following two sections we study pertur- 
bation dynamics adopting a simplified model where the 
photon-baryon plasma and cold dark matter are approx- 
imated by coupled fluids. The photon gas is assumed 



to behave as a locally isotropic fluid characterized by 
its density and velocity at every point in space and its 
pressure equals one third of the photon energy density. 
This is a good approximation before recombination when 
photons intensively scatter against free electrons and, be- 
cause of the Coulomb interaction between the electrons 
and baryons, the velocity of baryons is locked to equal 
the mean velocity of the photon fluid. We will arrive 
at these results consistently from the Boltzmann equa- 
tion and will consider the leading corrections to the fluid 
approximation in Sees. |v|-[v^. The effects of global cur- 
vature and of the cosmological constant should be very 
small on the scale of the acoustic dynamics before recom- 
bination and we do not include them in the discussion. 

Neutrinos contribute a large fraction (about 40%) of 
the radiation energy and require specification of their 
full phase space distribution even before recombination. 
We do not include a full treatment of neutrino pertur- 
bations in the present paper for the following reasons: 
The fluid model, corrected for photon diffusion and for 
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the neutrino contribution to the background energy den- 
sity as shown below, can describe perturbations in the 
early universe at least up to 5% accuracy, Ref. [||, when 
compared in its predictions of CMB anisotropy with full 
numerical calculations. The fluid description is adequate 
for baryons, dark matter, and photons before recombina- 
tion. Its simplicity and intuitive appeal make the fluid 
approximation an attractive starting point for applying 
perturbation theory in position space. The position space 
approach can also be used for accurate description of neu- 
trino dynamics and offers substantial advantages over the 
traditional Fourier decomposition, but that requires new 
constructions different from the spirit of the fluid descrip- 
tion, Ref. 1^. We postpone the phase space analysis of 
neutrino perturbations to a later paper. 

To account for the substantial contribution of neutri- 
nos to the background radiation density while excluding 
them cleanly from the perturbation dynamics, we con- 
sider a flctitious universe filled with photons, baryons, 
and cold dark matter (CDM), but without neutrinos. In 
this model universe the photon background energy den- 
sity, p^^°'^'^^\ equals the actual (physical) energy density 
of the combined relativistic backgrounds of photons and 
neutrinos: 

p(-°dcl) ^ ^(Phys) ^ ^(phys) ^ ^(phys) ^ ^ ^(phys) 

where 



Pi + Pv 



4/3' 



H -1 



is the fraction of the radiation energy density in neutri- 
nos that equals R^, ~ 0.408 for Nf ~ 3.04, e.g. from 
Ref. [||. 

The sound speed in the photon-baryon fluid is deter- 
mined by the ratio of baryon and photon energy densities 
as 



diPj + Ph) 



1 



adiab 



3 [1 + (3pb)/(4p^)] 



(7) 



The sound speed controls the dynamics of acoustic per- 
turbations via the length scale J Cgdr. To preserve this 
scale in our model, where is replaced by prei, we in- 
crease the baryon density by a factor {1 + R,^) over its 
physical value: 



Pb 



(model) 



= il + R.)p)^ 



(phys) 



(8) 



Finally, the total mean density of non-relativistic mat- 
ter pm, including baryons and CDM, is taken to equal its 
physical value: 



(model) 



(phys) 



(9) 



For eqs. (|2~|9|) to hold, the mean density of CDM in 
our model must be reduced slightly compared with its 



physical value: pc 



(model) 



„(phys) 

Pc 



R.pi"'''"'^. With 



these definitions, we arrive at a self-consistent two fluid 
model that preserves the important time scale of the 
radiation-matter density equality as well as the acoustic 
length scale J Cgdr. From now on, we drop the super- 
script "(model)". 



B. Dynamical equations 

The dynamics of perturbations in the fluid model is 
governed by the linearized Einstein and fluid equations. 
We give these equations in the conformal Newtonian 
gauge, Refs. |, |, and then der ive an equivalent but 
a more intuitive and easier to solve system of equations. 

Metric perturbations are induced by perturbations in 
the energy-momentum tensor T'V — X^a "^ai^ where a 
runs over all radiation or matter species in our model 
(a 7, 6, c). For every species, T^^ can be parameter- 
ized by an energy density enhancement 5a and a velocity 
potential Ua, denoted by W in Ref. [n4[. 



rpO 



~ (Pa+Spa) , 
{Pa + Pa) Vai , 
S} {Pa +Spa) ■ 



5pa = PaSa , (10a) 
Vai = -ViUa , (10b) 

(lOc) 



The stress T*j is isotropic for perfect fluids and the pres- 
sure perturbations equal Spj — i p^S^ and Sp^ — 6pi, — 0. 
Tight coupling between photons and baryons implies the 
equality of mean local velocities of the radiation and the 
baryon fluids. Sec. VI. Assuming adiabatic initial condi- 
tions, 



and 



Uh — Uy 



(11) 



before recombination. 

The linearized Einstein equations for metri c pe rturba- 
tions (f) and i/j in eq. (^ are given by eqs. (Al) of the 
Appendix. The last of eqs. (Al) states that for perfect 
fluids, when anisotropic stress is negligible, ip = <j). Then 
the remaining equations are 



.-H3-. 
a 



2« 
a 



%\^A^Ga'Y,^Pa. (12a) 

' a 

-(j)^AT:Ga^Y^{pa+Pa)Ua, (12b) 

a 

) = 47rG'a2 ^ (5pa , (12c) 



where "V" and "'" denote partial derivatives with respect 
to the comoving spatial coordinate r and the conformal 
time r. 

The fluid equations for density and velocity evolution 
in the conformal Newtonian gauge were derived in Ref. |^ 
and follow from the general formulae of Sec. |^. In the 
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tight coupling limit and our notations they are 



Sr. = 



6r 



■ Ur 



a 



1 

1 + (3p,)/(4p^ 



u, 

4p-y a 



(13a) 
(13b) 

(13c) 

(t> . (13d) 



The scalar gravitational potential (j) does not present an 
independent dynamical variable in addition to the densi- 
ties and velocity potentials because on any given hyper- 
surface of constant time the gravitational potential can 
be determined from the energy and momentum density 
perturbations on the same hypersurface via the general- 
ized Poisson equation 



47rGa 



5 Pa 



3— {pa +Pa)Ua 



following from eqs. (12a- 12b) 
We define 



Pm 
Prcl 



Pb 
Pm 



(14) 



(15) 



where Oeq — 1/(2.40 x 10"* ftmh^) is the scale factor value 
at the time of matter-radiation density equality and 
Pici = p-y hi our model. The ratio (3 is time indepen- 
dent. An important scale in our problem is a character- 
istic conformal time of the transition from radiation to 
matter domination 



(16) 



For the following ACDM model set of cosmological pa- 
rameters: i7i„ = 0.35, r^A = 0.65, Vlifh^ = 0.02, and 
h = i7o/(100kms~^ Mpc^^) — 0.65, the numerical value 

of Te is CTe ~ 130 MpC. 

With the above notations, the factor AnGa? on the 
right hand side of the Einstein equations becomes: 



2ri 



1 

ypn 



(17) 



If the dark energy is neglected in the early epoch, the 
Friedmann equation yields: 



1 



y 



y[T) 



T 



(18) 



The equality of matter and radiation densities, at which 

2/(Toq) = l, occurs at Tcq = 2(\/2- l)Te~0.83re. 

Dynamics of perturbations in our model depends on 
two characteristic speeds, given by 



1 



and 



dp 

d~p 



3 1 



adiab 



f3y) 



3 1 



(19) 



(20) 



P = J2aPa and p = J2aPa^ ^ee note [|5). Eq. m|) fol- 
lows from the sound speed definition of eq. ([^^since 

^(Phys)/^a^hys) ^ ^(modcl)^^(modcl) ^ ^^^^ ^^^^^^ 

speed, Cw, is not a true sound speed. It relates in- 
finitesimal pressure and density changes for perturba- 
tions of constant radiation entropy per unit mass of non- 
relativistic matter, ry, defined by 



77 = (5 ( In 



(21) 



= {1-P) ( -Sr-Sr 



We describe 77 by a dimensionless entropy potential a 
such that 



2t2 



(22) 



Using 0, 0, (7, a for independent dynamical variables. 



^ + (3 + 3cl 



3c2) 



3< 
Jrjy 



^a^yict-ci:^ 



(23) 



The first equation follo ws from the substitution of the 
left hand side of eqs. (12a, ^^) into the first line of 
eq. ( ^ and usi ng eqs. ( 
from eqs. (p^,[2l||l3|,[l2a| 
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12b) 



181). The second is derived 



Next, we replace the potentials (f> and cr by a pair 
of their linear combinations (j)r oc (f> + a/y and (j)c oc 
(cg/c^ — 1) <f) — a/y, which are chosen to diagonahze the 
second derivative terms in the system (|2^): 



i'c + Y.i=r Mci4>i + bci(f'i) = . 



The new variables 0^ and 
are normalized so that 



Then 



^-1 



and 



c2 


4> + 


(J 






y 












) 







(24) 

are uniquely defined if they 

0e . (25) 



(26) 

(27a) 
(27b) 



The matrices and bij are obtained by straightforward 
and somewhat tedious substitution of eqs. (If-H) in the 
system of eqs. (p3|). The result is 
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3 (c 



(i>c - 5 {cl - cl) - + 3- > 



4r2y 
3^ . , iVcl) 



3T2y2 



(2/c2)^4 + 3cg 

3t2?/2 

- 3 

'c = 



(28a) 
(28b) 



where = (1-/3) [2 + 2cl + 12^ + iclcl,{2 - icl - 3c^)] . 

Eqs. (^8|) are our principal equations of perturbation 
evolution in the model of two perfect fluids. We would 
like to stress that these are causal hyperbolic (wave) 
equations in contrast to "action at a distance" type el- 
liptic equation ([T^). This difference from the original 
Einstein equations enables straightforward numerical in- 
tegrati on of eqs. ( p8| ) directly in real space. 

Eq. ( p8a| ) shows that perturbations in the radiation- 
baryon plasma, described by (j)r, propagate as acoustic 
waves with the sound speed Cs- As the radiation wave 
passes by, its gravity perturbs locally the c old d ark mat- 
ter as described by (jfc and (j)c terms in eq. (28b), and the 
CDM potential 0c starts evolving. 



described in terms of 



Then by eqs. (k5 



e — S + 3 — u 

y 



e^Mv2(y0) 



(32) 



(33) 



and e satisfies the conservation equation e + V • = 
with = — Vu. 

The linear combination on the right hand side of 
eq. ( p^ ) remains invariant under redefinition of constant- 
time hypersurfaces, f = t + a(r,T), and spatial coordi- 
nates, f = r + V/3(r, t), when in the linear order 



C. Density and velocity fields 

In order to analyze the CMB temperature anisotropy 
or matter distribution in the universe, the potentials (j)r 
and (f>c should be related to perturbations in density and 
velocity fields. The corresponding equations are derived 
in this subsection. 

It is convenient to replace and by the following 
variables linear in y: 



Oa= Oa a 

Pa 



(34) 



(For any a^O or f3 ^0, the metric ds'^ will no longer be 
conformally isotropic in space as it is in eq. (|6|).) Using 
the equation of energy conservation. 



Pa = -3 - (pa +Pa) , 

a 



we can see that 



A: 



1 , 3 



For the quantities 

.E^Pa 



J2iPa+Pa)Ua 



(29) 



(30) 



the entropy perturbation 77 of eq. (|2l|), and for an addi- 
tional variable 9 defined below, from eqs. (p^ p^,[T7|) we 
find 



■5r + {l~P) 5, 



2< 
3 



iB 
3v 



Ur + {I - (3) Uc 



3f (2/0)- 



77=(l-/3)(|<5,-<5c) =^V2a 



(31) 



Since by the last two equations + V • jr, = with = 
— V0, physically, is a potential for the entropy current 
density. 

The above equations become particularly symmetric 
in and a if the perturbation of energy density is 



Hence the invariant variable e of eq. (|3|) equals the en- 
ergy density perturbation 6 in the gauge where u is iden- 
tically zero, i.e. the combined momentum density of 
both fluids vanishes. Such a gauge is a generalization of 
the comoving gauge to our multicomponent system, and 
eq. ( ^2[ ) generalizes Bardeen's gauge invariant variable e 
introduced in Re f. |2| . 

Using eqs . (^-^2D and the definitions of (jjr and 0c hi 
eqs. (p5|-p6|), it is straightforward to find: 



3y . 



Sc 



3 3y , 

4 A 



where 



(^ + ^) = ¥ivMI#0 , 

^ r;=M_l_v2(y0c), 



(35a) 



(35b) 



and 



- 3y 



A(l-/3) 



3y_ 



B 

^(1-/3) 



2< _1_ 

3 1-/3 



(#c)' 



(A/By 

(A/B) 



(35c) 
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In addition to determining the density and velocity 
perturbations, eqs. (^5| ) also help to establish the physical 
meaning of the potentials 0^ and (j)c- From eqs. ( 35b ), 
remembering eq. (17), 



Thus, the potentials cjjr and (j)c are the metric pertur- 
bations generated individually by the (comoving gauge) 
disturbances in the radiation-baryon plasma and CDM. 

The following conformal Newtonian gauge quantities 
contribute to the primary CMB temperature anisotropy 
in the tight coupling regime, Refs. |16: |l^: j6j. + 4>, 
Ur, and 0. Of these, Ur is related to 0^ by the first 
of eqs. (35c) and (j) is given by eq. (|25|). We can compute 
j6r + (f> from Ur and </> as 



1 3 

-6r + (p = Ur + [yUr) - y<j) 



(37) 



This equation follows from eq. (13d) 



III. 



GREEN'S FUNCTIONS AND INITIAL 
CONDITIONS 

A. Defining equations 



The evolution of cosmological perturbations in the lin- 
ear regime may be studied in position space by super- 
imposing their individual Green's functions. For the 
Green's functions one can take any complete set of so- 
lutions of eqs. ( p8| ) satisfying suitable initial conditions. 
We consider two types of initial conditions: An initially 
point-like perturbation at the origin, e.g. 



</,(3)(r,r) 



as T 







(38) 



where ^^■'(r) is the product of spatial coordinate delta 
functions SY){x)Suiy)Su{z) , and a sheet-like perturbation 
that is produced on the whole (0, y, z) plane and is inde- 
pendent of the y and z coordinates, Ref. H], e.g. 



(/.W(a;,r)-.<5D(x) 



as T 







(39) 



The initial conditions in the form ( p^ ) will prove to be 
the most attractive for describing perturbation evolution. 

All the fields 4>r, ipc etc. are functions of three 
spatial coordinates x, y, and z, in addition to r. If, 
however, a configuration of fields has no initial depen- 
dence on y and z, as in eq. (39), the fields will remain 
independent of y and z for all time. For such a configu- 
ration one needs to solve the evolution equations (^8|) in 
the X and r coordinates only. Thus a Green's function 
for the evolution equations with the plane-parallel initial 
conditions, eq. (p^, is also the Green's function for the 



same equations considered in one spatial dimension. We 
denote these Green's functions by the superscript "(1)", 
and use the superscript "(3)" to distinguish the perturba- 
tions originating from a point disturbance, as in eq. (|3^). 

When setting the initial conditions, one should also 
specify the ratio of perturbation magnitude for differ- 
ent species, in our case the radiation and the dark mat- 
ter. In this paper we consider only adiabatic initial con- 
ditions, which are favored by the present experimental 
data, Refs. @ [l9[ 1^, Jl, Hi, and by the simplest infia- 
tionary models, Refs. [ESl: 



y<p 



as 







(40) 



(A Green's function approach to isocurvature initial con- 
ditions, arising naturally in topological defect models, 
Ref. Q, has been considered in Refs. ||^, |ll|.) 

Given the initial ratio of perturbations in various 
species, e.g. as implied by eq. (^0|), the subsequent evo- 
lution is completely described by a single set of Green's 
functions, either one-dimensional or three-dimensional, 
for the relevant variables, in our case (l)r and (pc. To il- 
lustrate the use of Green's functions, consider the total 
gravitational potential (p. The time dependence of its 
Fourier modes is given by 



0(k,T)=0(^)(fc,T)A(k) 



(41) 



Here, A(k.) is a primordial potential perturbation created 
by inflation, Ref. ||2^ , or another mechanism of perturba- 
tion generation, e.g. Refs. [^J, 26|, and 4'^'^\k,T) is the 
Fourier space transfer function 



0(^)(fc, 



or 



\k,T) = / (fxe-'^-''(l)^^\r,T) 



(42) 



(43) 



Comparing eq. (^) and eq. ( p2[ ) 
types of Green's functions as 



one can relate the two 



r+oo 
J\x\ 



(44) 



Eq. (^) says that in order to find the magnitude of the 
perturbation at a distance x from the plane of the initial 
excitation ( p9| ) , one should add up the perturbations orig- 
inating at all the points on the plane which are separated 
from the considered point by r = ■y/aJ^ + j/^ + z^. Con- 
versely, a three-dimensional Green's function can always 
be obtained from a plane-parallel one by differentiation: 



)(r,T) 



dr 



(45) 
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B. Radiation era solutions 



When radiation dominates the background density, 



1/3 and c2 



y ~ - < 1 



(l/4)(l-/3)y.. 



For the adiabatic initial conditions ([40D, eqs. ( |27D 
(pc'^yipr during the radiation era. Retaining in eqs. 
the leading in y terms only, 



nve 

m 



0, = 3(1-/3)2/ 



1 

17' 



Their non-singular one-dimensional Green's function so- 
lutions satisfying initial conditions (Is^) and (Eoh are 



4 [CsTy 



'(x,t) 



9(1-/3)2; 



3(c,r)2 



2c,t|2:| — X 



(c.r)2 |x| 



4(c.t)3 

c,r -Ixl) , (47b) 



where 9 is the Heaviside step function: 9{x') — 1 for x' > 
and 6'(x') — otherwise. The Fourier transforms of 
these Green's functions are the "growing mode" transfer 
functions 



4^\k,T) = 3 



sm z cos z 



9(1-/3)2; 



z = kcsT , (48) 



sin(z) — z In z + C — ci(z) 



2z3 



where C = 0.5772... is the Euler constant and 

cos z' — 1 



■dz' ^\iiz + C + 



dz' 



ci(z) 



is the cosine integral. 



IV. PROPAGATION OF PERTURBATIONS 



The one dimensional Green's functions for eqs. (|2£ 
at a later time can be obtained by numerical integration 
of these equations starting from the radiation era solu- 
tions (^^ at some rinit <C t^. The applied numerical 
methods are described an d com pared with Fourier space 
calculations in subsection IV C| . The results of the inte- 
gration up to the time of hydrogen recombination at red- 
shift z ~ 1.07 X 10'^ are presented in Fig. |l|. The original 
delta function perturbation has separated into left-going 
and right-going waves, whose evolution spreads ^r^"* over 
space and diminishes the magnitude of the radiation dis- 
turbance. The gravitational potential hill of the radiation 
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FIG. 1: Green's functions for the potentials (solid) and 
(/>c (dashed) at the time of recombination Tree — 2.15 Te with 
an = 0.35, = 0.65, a/i^ = 0.02, and h = 0.65. 
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FIG. 2; Time evolution of the Green's functions of Fig. ^ 
describing the gravitational potentials for radiation (left) and 
CDM (right). The plots, from top to bottom, correspond to 
the conformal time values 0.1 re, Te, and Tree — 2.15 r^. 



causes outward-directed gravitational forces which expel 
the CDM away from a; = 0. Three snapshots of the time 
evolution of ^r"'^' and (\)^c^ are shown in Fig. ||. Only the 
range a; > needs to be calculated since the potential 
Green's functions are even functions of x. 

The Green's functions are identically zero beyond the 
acoustic horizon, \x\> S{t), where 



Sir) 




(49) 



{S{t) ~ T I \f?> when t ^ Tg). In a more careful treat- 
ment one will find weak precursors extending beyond the 
Green's functions acoustic fronts up to the particle hori- 
zon X — ±CT. The precursors arise because of partial 
photon, and in a full treatment neutrino, free streaming 
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FIG. 3: Green's functions of the radiation density enhance- 
ment (top) and the CDM density enhancement 5c (bottom) 
at the time of recombination. The vertical spikes represent 
the delta- function singularities of irix) at its wavefronts and 
of (5c (a;) at the origin. 



with the speed of light. Postponing detailed discussion of 
the photon diffusion until Sec. Vl and of the neutrino free 
streaming until a separate paper, we continue the analy- 
sis of perturbation dynamics in the perfect fluid model. 

The density and velocity perturbations of the photon- 
baryon and CDM fluids are determined from the (\)r and 
(^c potentials by eqs. (|35|). The corresponding Green's 
functions are plotted in Figs. ^ and |[ Characteristic 
features in the density and velocity Green's functions 
and their connection to co smological parameters are dis- 
cussed in subsection IV D The delta function spikes at 
the acoustic wavefronts of and Wr^"* in the fluid model 
are fully calculable analytically and are considered in sub- 
section IV B. 



A. Sum rules 

Sum rules are simple, easy to apply integral relations 
that offer a powerful tool for checking analytical formu- 
lae or numerical calculations. The sum rules are more 
general than the fluid approximation: they hold in the 
linear regime regardless of how complicated the internal 
dynamics of perturbations may be. They offer a simple 
method for estimating the accuracy of numerical calcu- 
lations of the Green's functions. In addition, the sum 
rules give insight into the connection between the posi- 
tion and Fourier space descriptions. They become in- 
dispensable for setting the initial conditions of Green's 
functions when one considers perturbations to the full 
phase space distributions, Ref. 0. 



The idea of the sum rules is to use eq. (^2|) with fc = 0: 



0(T)(/!c = O,r) 



da;0(^)(a;,r) 



(50) 



The time evolution of the k — ^ Fourier modes, the mag- 
nitude of which wc denote by upper case letters, e.g. 
0'^-^-' (fc =0, r) = <i>(T), satisfies easily solvable ordinary 
differential equations. The fc-space equations at fc = are 
especially simple because, given the adiabatic initial con- 
ditions, all the parameters characterizing relative motion 
of different matter and radiation components are unper- 
turbed at all time on infinitely large scales. For example, 
in our case of the radiation-CDM model. 



(51) 



Trivially integrating first of eqs. ( |2^ ) over x from — cxd 
to cxD, we obtain the equation for the gravitational po- 
tential in the /c = mode: 



y 



(52) 



With given by eq. (20), the linear independent solu- 
tions of this equation are 

9y 9y2 9y3 ' 



Taking the linear combination of /i and /2 satisfying the 
initial condition <I>(t — > 0) = 1, we obtain the sum rule 
for the gravitational potential: 

$(r) ^Jdx <t>'^'^ (x, r) = ^h+^^h. (53) 

The sum rules for (pr and (pc then follow immediately 
from eqs.(p7|,p|,Fll): 



(54a) 
(54b) 



The sum rules for the radiation velocity potential Ur 
and the combination j5r + (p are relevant for CMB tem- 
perature anisotropy. The first sum rule sim ply follows 
from the first of eqs. (35c) and from eq. (54a): 



C/.EE j dxu(^\x,T) - ^^Ml^y,^), (55) 

The second can be derived starting from eq. (^7|) and 
applying formulas of this subsection, including eq. (^2|), 

iA. + <i> ^ j dx(^^e)+<t>(')^ = 

^U.^-^M^(^-A' ■ (56) 
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FIG. 4: Green's functions of the velocities Vr 



■ Vwr.c (left) and of the velocity potentials tt^.c (right) at the time of 



recombination. Solid and dashed lines are for the radiation and CDM respectively. 



Finally, we give the sum rules for the radiation and 
CDM energy density enhancements. They follow from 
eqs. (35a -35b) and the second of eqs. (|3l|): 



A. = -A, = --|i^ (y$) 



(57) 



with = J dx (5r^'' and = J dx S, 



(1) 



B. Wavefront singularities 

We can obtain exact analytic formulae describing 
Green's functions at the acoustic wavefronts |a;| — S{t) 
in the fluid model. As one will see below, the wave- 
front spikes in the density and velocity Green's functions 
play a significant role in the acoustic dynamics. For the 
wavefront analysis, it is convenient to factor the radiation 
potential as 

Mx,r) = C{T)d{S{T)~\x\,T) (58) 

where for the perfect fluids, with x' = S{t) — \x\, 

d(x'<0,T) = 0, Ad(+o,r) = l, (59) 

i.e. d{x' ,t) ~ x' 9{x') as x' — ^ 0. Substituting these 
equations in eq. (28a), setting |a:| = S[t), and taking 
into account that at the wave fronts 



we find the following simple relation for C(r): 



2(7- 



C = 



(60) 



(61) 



Here we have applied the formula 
3c^^ = ^ + 2^ 

y y cs 

for the sound speed given by eq. (|l^) . Integrating eq. ( |6l| ) 
and norm alizi ng the result to agree with the radiation era 
solution (47a) in the y limit, we obtain: 



2r2y2 (3c2) 



In particular, at a; = ±S'(r), 



3\/3 



2r2y2 (3c2)i/4 



(63) 



The values of velocity potentials at the wave fronts 
follow from eqs. (p5q): 



,(1) 



3\/3 



,3/4 



,(1) 







(64) 



Because 4>r and Ur are identically zero when > S{t), 
the velocity Vr = — Vt tr and the radiation density 5r, 
given by eqs. ( p5a -35b), both contain a delta function 
singularity at the acoustic horizon: 



'5!.:U = ^!^!U = 3(3c?)'/'fe(|:r|-^) 



(65a) 



^(3c2)'/%ign(x)fc(N-5) . (65b) 



The dynamics of the shock-like singularities in the 
Green's functions is linear, even if nonlinear effects would 
take place in the propagation of a real shock wave in 
the photon-baryon plasma. Indeed, the actual cosmo- 
logical perturbations are the convolutions of the Green's 
functions with the smooth primordial potential field. As 
long as the overall perturbations are small, they remain 
as regular as the primordial perturbation field and their 
dynamics is well in the linear regime. (On very small 
scales, there may exist physically relevant nonlinear ef- 
fects, Ref. ||2^, unrelated to the above Green's function 
singularities. These effects are absent when photons and 
baryons are tightly coupled.) 



C. Numerical integration 

From the hyperbolic (wave) nature of eqs. (p8|), the 
values of the potentials (pr and 0c at a point (x, t -|- At) 
are uniquely determined by their past values within the 
interval [x — CsAt, x -I- CsAr] at time t. Cusps at the 
potential wavefronts will not be distorted by numerical 
errors during the evolution if discretization intervals in 
space and time are related so that c^Ar is a multiple of 
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Aa;. In particular, for c^Ar = Ax, the second derivative 
terms in eqs. (p^ ) may be approximated by the second 
order scheme 



these calculations as estimated by the sum rules is 0.19% 
and 0.05% respectively. 



~ '^[(j)r{x,T + At) + (j)r{x,T - At) 



Mr) 



(At) 



(At) 



^{x + CsAt, t) - (f>r{x- CsAt, t)] , (66) 0, 



[0c(r + At) + 0e(r - At) - 2^,{t)] 



In practice, we use the horizon radius S{t) for the time 
evolution coordinate and take AS = Ax. 

The first time derivatives of c/fr and (j)c can be approx- 
imated by an implicit second order scheme 



(r) 



1 

A^ 



[0(t + At) - 0(t)] 



At 



(67) 



An explicit scheme should be avoided as numerically un- 
stable. With the substitutions (|66H67|), eqs. ( p8| ) become 
finite difference equations that can be straightforwardly 
solved for the values of (j)r and 0c at {x,t + At). The 
resulting finite difference scheme is of the second order 
and it preserves the form of cusps or discontinuities in 
propagating waves. 

We start the integration from the radiation era solu- 



tions (47b) at some small initial time, e.g. Tinit — 10~ Tg. 
To follow the evolution over a large dynamic range, the 
spacings of the space and time grids are increased by a 
factor of 2 every time the extent of the acoustic hori- 
zon S doubles. Then the number of discretization points 
N = S/ Ax used to represent the Green's functions at 
all later times remains in the range iV,„in < N < 2Njnin 
where A^min is the initial number of points within the hori- 
zon. The total number of c/fr and (/)c evaluations needed 
to evolve the potentials to a final time Tfin is of the order 

of iV2 ln(Tfi„/Tinit). 

Whenever a Fourier space transfer function is needed, 
it can be obtained from the corresponding Green's func- 
tion using the Fast Fourier Transform (FFT) algorithm. 
The FFT can be efficiently applied to Fourier integrals 
as described in Ref. ||2^ with its CPU time scaling as 
A^log2 N. Due to the finite extent of the Green's func- 
tions, the related transfer functions oscillate in k and t 
requiring in general a larger number of evaluations for 
their direct Fourier space calculation with the same ac- 
curacy. 

prove to 



The sum rules introduced in Subsection IV A 



be a highly efficient debugging tool available for calcula- 
tions in position space. The sum rule relations can also 
be used to estimate the numerical error of computations. 

Calculations of the Green's functions in Figs. |l| and |^ 
proceeded from Tjnit — 10~^Te to Tree — 2.15 Te. The 
corresponding CPU time for our Pentium IV PC is 0.04 s 
with the minimal number of points in the spatial grid 



D. Discussion of Green's function features 

The plane parallel Green's functions for the potentials 
and (j)c, plotted in Figs. Q and |^, are continuous func- 
tions of X and t. In the fluid approximation, the radiation 
potential is characterized by discontinuous first deriva- 
tives (kinks) at its acoustic wavefronts as described by 
eq. (p3). The CDM potential has a central cusp refiect- 
ing the initial repulsive singularity in the gravitational 
potential i^f^''; this cusp is preserved because the CDM 
particles have no thermal motion. In fact, once the pho- 
ton energy density becomes negligible in the matter dom- 
inated era, the potential (pdx^T) stops evolving, as may 
be seen from eqs. ( |2^ ) in the limit y ^ 1. 

The discontinuities in the potential derivatives give 
rise to delta function singularities in the corresponding 
density and velocity transfer functions, following from 
eqs. (p5| ) and plotted for the time of recombination in 
Figs. I^and ^. The singularities are visualized in the fig- 
ures by the vertical spikes aX x = ±S for 5^ and Vr and 
at X = for 5c. 

The wavefront singularities in 5^r^ {x) and Wr^"* {x) are 
described analytically by eqs. (^5|). Their significance in 
perturbation dynamics may be characterized by the delta 
function contributions to the sum rule integrals for the 
corresponding Green's functions. For the radiation en- 
ergy density, we find that in our ACDM model at recom- 
bination / dxSi^Ki^^/ J dxSi^'^ ~ 2.7/(-2.3). Thus the 
delta function singularities play a major role in the per- 
turbation dynamics even on the largest scales. Their rel- 
ative contribution to the CMB anisotropy even increases 
on smaller scales because the oscillation amplitude of the 
delta function Fourier transform does not fall off with k, 
as opposed to the oscillations in the non-singular term 
Fourier transforms. We have shown in Ref. ^ that in 
real space the singularities give rise to a localized feature 
in the CMB angular correlation function C{d). 

The delta function singularities acquire a finite width 
when photon diffusion, or Silk damping, Ref. p9[ , is in- 
cluded into the analysis. When the photon mean free 
path is small compared with the scales considered. Silk 
damping can be approximated as 



<Sf)(fc) 



(T) 

r fluid 



Refs. [BO Pl p2 . As foimd in Ref. [02 



4W - g 



14 , , 

15 (^^^ 



(3, 



,2 V 



cdT 



(68) 



(69) 



64 and 0.1s with N„ 



128. The accuracy of 



where Ue is the proper number density of electrons and 
aT is Thomson scattering cross section. Silk damping 
broadens a delta function singularity in a position space 
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FIG. 5: The position space transfer function for the intrinsic 
ijS^) plus gravitational redshift contributions to the CMB 
temperature anisotropy AT/T at recombination assuming 
tight coupling between photons and baryons. The values of 
cosmological parameters are taken as fim = 0.35, SIa = 0.65, 
and h — 0.65. Baryon drag effects are evident from comparing 
the solid, Qbh^ = 0.02, and the dashed, Q,b = 0, lines. 



Green's function to a Gaussian of finite width 

\2' 



5t>{\x\^S) 



1 



exp 



4x| 



(70) 



as obtained by applying the filtering (68) to the con- 
stant Fourier transform of a delta function. We red- 
erive eqs. (|68|-|69|) within our position space formalism 
in Sec. 0. 

The CDM density perturbations, with the Green's 
function shown in the bottom panel of Fig. 3, eventu- 
ally seed the formation of galaxies. The central spike in 
the figure arises from V^(/<c in the second of eqs. (35b). 
It is negative because of the repulsive sign of the initial 
gravitational potential peak, as shown in the top row of 
Fig. |. The spike is surrounded by positive tails because 
the CDM pushed out from x = piles up into the region 
between a; = and the acoustic wavefront. 

To understand the physics of the dip in (x) at a; = 
in Fig. ^, let us examine the intrinsic radiation temper- 
ature perturbation (5^/4 corrected for the temperature 
redshift by the local metric perturbation at x: 



■ 5r + ' 



(71) 



The plane-parallel Green's function for Aes at recombi- 



(1) 



at 



nation is shown in Fig. ^. The smoothness of A^^^ 
x = in the absence of baryons (dashed line in Fig. ||) 
can be understood as the condition of thermodynamic, 
Ref. | |3^ , and hydrostatic photon equilibrium near x = 0. 
(Obviously the equilibrium does not apply at the acoustic 
wavefronts, but it does hold as a; ^ where the elapsed 
time is many sound-crossing times.) Linearizing the rel- 
ativistic equation of hydrostatic equilibrium, 

V,p+(p + p)V,0 = O , (72) 

for the case of the radiation fiuid containing no baryons 



V,; I - (5^ 



4 
3 ' 



= - V.Aeff = 



(73) 



Thus, if ph = 0, Acff has zero gradient in equilibrium. 
If ilfc 7^ 0, on the other hand, the hydrostatic eq. ( |7^ ) 
applied to tightly coupled radiation-baryon fluid gives 
ViAcff = -(3/4)/3yV»(^. The positive cusp of the CDM 
gravitational potential becomes a dip of CMB tempera- 
ture. The dip appears because the heavy pressureless 
baryons are repelled by the cusp in the CDM gravi- 
tational potential (^^ from the a; = region and they 
drag the coupled photon fluid away from the origin. In 
eq. (p7|), rewritten as 



'■Mc. (74) 



only the last, 0^, term on the right hand side has a 
cusp at X = 0. The contribution of this term and thus 
the magnitude of the dip in Aoff are roughly propor- 
tional to (3y{T^cc) oc ilbh^. Another effect of the baryonic 
constituent on the radiation perturbations evident from 
Fig. ^ is the decrease of the sound speed and thereby the 
reduction of the acoustic distance compared with a pure 
photon gas. 



V. PERTURBATION DYNAMICS BEYOND 
THE FLUID MODEL 

Now we will distinguish between photon and baryon 
perturbations and will not assume the fluid approxima- 
tion for photons. We will continue to exclude neutrino 



perturbations from our model as described in Sec. [I A. 

The linear evolution equations for the photon and 
baryon velocity potentials follow from linearized Boltz- 
mann equation and are given in the Appen dix a s the 
second of eqs. (A12) and the second of eqs. (A15). In 
these equations, the terms proportional to the inverse of 
the mean conformal time of a photon free flight, 

T^^ = aUeCTT , (75) 

are due to Thomson scattering between photons and 
baryons. The Thomson scattering damps the velocity 
potential difference, for which we define 



Ud = Ub 



(76) 



However, the scattering does not affect the overall mo- 
mentum density of the plasma. Defining a momentum- 
averaged velocity potential of photons and baryons 



Ur 



{p.y+P't)+Pb B 



3py 



■ Ub 



(77) 



where B is g iven by eq. (p9[), from the velocity equations 
in eqs. (|A12 , |A15|) we find: 



B \4"^ 



Ud 



Mi 



7 



J_ 4B_ 

Ta 3f3y 



(78) 
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Substituting = Ur — Ud and Uh = 4 Ud + in 
the density evolution equations in eqs. ( A12 , A15| ), and 
defining 



5d 



(79) 



to photons: 



2r2 1^ T3 



2Ti 



(83) 



we obtain 



6d ^ V^Ud . 



(80) 



The cold dark matter density and velocity perturba- 
tions evolve as 



Sip 



(81) 



The conformal Newtonian gauge potentials (/> and ip, de- 
fined by eq. (^), can be determined non-dynam ical ly from 
the first, second, and fourth relations of eqs. (Al): 



1p -(/)-- 



2r2y 
6 



B 



y 



y 



y 



9 9 "7 
Til/ 



(82a) 
(82b) 



The other can be the entropy potential a from 
eqs. (E2lElh that now equals 




5-1 — 5c 



pSd 



(84) 



The time derivatives of V^ad or V^u equal a full V'^ 
of a certain linear combination of velocity potentials Ua- 
Differentiating both sides of eq. (p3) a nd eq. ( p4|) with 
respect to r, remembering eqs. (|80|- 8l|), and lifting V^, 
one obtains: 



3 



(1-/3) {Ur-U,)^^Ud 



(85) 



whe re fo rmula ( |l7] ) was applied. The variable u 
eq. (B2a) is defined as previously by eq. (BOh. 



Two independent variables are needed to describe the 
species number density perturbations jS^, 61,, and 6c rel- 
ative to each other. For one such variable we use a poten- 
tial generating the baryon density perturbation relatively 



Taking another time derivative of both sides of 
eqs. (p5| ), using the dynamical iia equations from 
eqs. (^H), and expressing all the resulting density 
and velocity perturbations in terms of 'S/'^ip, as given by 
eq. (|23), V^aa, eqs. (|8l||), and &a, eqs. iM), 
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(86a) 
(86b) 



Similarly to the fiuid case, these equations are supplemented by the formula following from the first line of eq. (pi]), 
eqs. (|2[|l^, and the Einstein equations (Al): 



(3cU2)^V'+^' 



4r2y 



7 



cr 2tt. 

V + - + -2— 2 2 

y ctjTey 



(86c) 



Here (j) is expressed via ip and ■n-y using eq. (82b). The 



set of eqs. (86, B2b) is not closed as long as tt^ remains 



an independent variable. 

We split the total potential ip into potentials ipa gen- 
erated by the individual species: 



■0 = V'7 + V^c + V'd 



(87) 



The dynamics of the potentials ipa may be described by 
coupled wave equations provided the decomposition (p7^) 



is performed according to the two following requirements: 
First, on small scales, when the velocity term in eq. (82a) 
is negligible, every V^V'a is proportional to the corre- 
sponding 5a and is independent of the density of the 
other species. Second, the potentials ipa are given by 
linear combinations of ip and the entropy potentials a 
and Gd- These linear combinations are unique and are 
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easily found to be 



ci \ V 



ci J \ciJ y y 



f3(Td 

y 



The corresponding Laplacians are 



2T2y V y 

3 



2T2y 
3 

2^ 



(1 - P)Sc 





3 - M 









(88b) 
(88c) 

(89a) 
, (89b) 
(89c) 



The evolution equations for the potentials ipa a-i'e ob- 
tained by substituting eq. (Is^) and 



y 



, crd = y^^d (90) 



in eqs. (B^, 82b). Upon the substitution, one finds a sys- 
tem of coupled second order wave equations generalizing 
eqs. (p^). For simplicity, here we write these equations 
omitting the terms with no or only first derivatives of 
the potentials ipa and tt^ unless a term contains the large 
damping parameter t~^. The neglected terms contribute 
only to the dynamics on large scales A ^ r, when the 
fluid eqs. (|2^) can be used. With the remaining terms, 
one obtains that on scales A <C r 



V'7 

4 



9/32/c2r, 



ipd 
■ipd 



Wycl ^ 

4 

. 



27r^ 



(91a) 

,(91b) 
(91c) 



In general, eqs. (|9^) should be completed by the evolution 
equations for the higher multipoles of the radiation phase 
space density f-yi>2 and polarization multipoles g-yi that 
are given in the Appendix. However, in the following 
section we find that in the next to the leading order in 
the photon-baryon coupling Tc, 7r-y = ^ /-y 2 is fully deter- 
mined by ipj and ip^. Then, neglecting special features of 
neutrino dynamics, the above equations provide a closed 
system. 



find the 0(tc) leading corrections to the fluid approxi- 
mation. 

First of all, iid equation in eqs. ( [T^ ) and all the equa- 
tions for f-fi>2 and g-^i in the Appendix have the generic 
form 



fi 



— {.fi - b) , 

rTr 



(92) 



where r is a positive number and a and b are some 
linear combinations of variables other than //. In the 
tight coupling regime, <^ r, the evolution given by 
eq. ( |9^ ) will quickly, over a time of order Tc, drive // to 
// ~ b + rTc{a — b), as is evident from the explicit solution 
of eq. (H) 



h = 5 + rdr'(a-5)exp(-/;^ 
= b + rT,{a - i}) + 0{t^) . 



(93) 



Thus at tight coupling the dynamical equations of the 
form ( |9^ ) reduce to algebraic equations /; ~ 6 up to 0{t^) 
or, if b in eq. (^) is absent, /; ~ rrcfl up to 0{tc). 

Applying this observation to the evolution of the pho- 
ton polarization averaged multipoles f^i in eqs. (A12) 
and to the polarization difference multipoles g-yi in 
eqs. ( A14| ), we see that each / > 3 multipole is suppressed 
relative to the corresponding lower multipole by an extra 



r^, we then find that 



power of Tc ■ Since 1 = ^u^ 
III ^ Tc"^ for I > 1: 57O ~ Tc, g^yi ~ T^, and g-^i ~ 
for Z > 2. The equations for /^2, 5^0 ^^'^ 57 2 reduce to 
the following algebraic relations: 



A 2 — (A 2 



2 

- 57O + .972) +Tc-fjl 



^7 — 2 ( A 2 + ^7 + 57 2) , 
572 — (/72 + 37O + 372) 



(94) 



where the omitted terms are 0{t^). Similarly, from the 
second of eqs. (|78|), 



Ud 



Wy 

AB 



1 



Resolving the linear algebraic system (£^ 
and then remembering eq. (95) we find: 



(95) 

|) in terms of /-y 1 



1 



72 



4 , 
rc-/7i 



16 

Tr U, 

45 



■7 



16 

Tr Ur 

45 



(96) 



Substitution of results (^,^) into the right hand side 
of 5y and iir equations in eqs. (3^,|7^) gives 



VI. TIGHT COUPLING: NEXT TO THE 
LEADING ORDER 

Now we can consider systematically how the general 
formulae of the previous section reduce to the fluid equa- 
tions in the limit of tight photon-baryon coupling and 



1 



( Mv\ V 

\iB J y 



Comparing these equations with their perfect fluid ana- 
logues, eqs. (13c -13d), one can identify two qualitatively 



15 



new terms appearing in the 0(tc) order. The V^^-,. term where 
in the first equation corresponds phenomenologically to 
heat conduction and V^itr. in the second to bulk viscos- 



ity of adiabatic scalar perturbations in the photon-baryon 
plasma. 

If Thomson scattering did not partially polarize the 
scattered radiation, in place of eqs. (Q) we would find 
fj2 — I 1 . This is 25% less than the result ( |9^ ) 
and would wrongly yield a 25% smaller value of the ra- 
diation bulk viscosity. The relevance of fluctuations in 
photon polarization for dissipation of the perturbations 
in photon-baryon fluid was pointed out in Ref. [3^ . 

The remaining undetermined quantity describing the 
relative motion of photons and baryons is their num- 
ber d ensi ty difference S^, which enters the Poisson equa- 
tion (B2a). We will calculate the related potential ad de- 
fined by eq. (|8^). Within the 0{tc) accuracy considered 
here one can substitute the perfect fluid expressions ( |35| ) 
into the right hand side of eq. (p5|). Then replacing (j)r 
and 4>c by (j) and a and using the second of eqs. (^3|) to 
eliminate {(j) + a/y), we find 



2r3 



Ud 



3/3 



3 4(1-/3) 



(97) 



The relation (Jd — 3/(2Tg)urf from eq. ( |85| ) can then be 
integrated from r = 0, with cr^l^^Q = for adiabatic 
perturbations, to a given t: 



Od 



3/3 



4(1 - /3) 



(98) 



Thus to first order in Tc one can continue to de- 
scribe the perturbation dynamics by only t wo in depen- 



dent fields 



and 



■\\)c- By eqs. 




2b|), these 



fields reduce to our earlier definitions ( |27| ) in Tc ^ 
limit. In the 0(tc) order, the potential il^d is related to 
0r and 0c via eq. ( ^8c] ) and eq. (|9^), in which a may be 
substituted by the 0(t°) order expression (p6|). 

Since ti-^ in eq. (|9q ) is an 0(tc) quantity, Ur on the right 
hand side of eq. ( |96|) may also be expressed in terms of 0,. 
using the leading order relations (|35cD . When the corre- 
sponding expression for tt^ and eq. (|98|) are used to close 
the dynamical equations of the previous subsection, one 
finds two coupled differential equations for 0^ and 0c that 
contain all the terms of eqs. (2^) and additional terms 
proportional to Tc. (Note that the third time derivative 
of a arising from '&d may be reduced within our 0(rc) 
accuracy by relations (^ to terms containing a lesser 
number of time derivatives.) When Tc ^ r, the 0(tc) 
terms containing three derivatives of the potential may 
still contribute appreciably to small scale variations of 
the potentials when their higher r or spatial derivatives 
become increasingly large. Retaining only the second 
derivative terms with O(t^) c oefficien ts and all the third 
derivative terms, from eqs. (|91a, 91c) we arrive at the 



following system, valid on small scales A <C t: 



2Tc.g V^0r , 



(99a) 
(99b) 



1 

52 



(100) 



The last term in eq. (99a) and so the expression ( |100| ) 
rece ive contributions from both t\}d and -n^ terms in 
eq. (|l^. 



The V^0r term in eq. (99a) describes the famous Silk 
damping of perturbations in the photon-baryon plasma 
on small scales, Ref. [|9|. In momentum space, the dis- 
persion relation imposed by eq. (99a) on a plane wave 
0r = exp (ik • r — Ilot) is 

+ 2iTcgk^uj - fc^cf = . (101) 

The solutions to first order in Tc are uj — ± kc^ — «7 with 
the damping rate 

7 - Tcgk^ ■ (102) 

This rate coincides with the result of Ref. |^^, where 
derivations were done with a different approach that also 
took into account the polarizing property of Thomson 
scattering. 

One can quantitatively describe photon diffusion at the 

sharp wavefront of 4'r^\x,T), as it evolves from Tinit — > 
to a certain finite r, by considering the full eq. (28a) with 
the third derivative term 2rcgV^0r added to its right 
hand side. We look for a solution in the wavefront region 
|S'(r) — |a;|| ^ T using the ansatz (|5^): 

0(1) = C(r) d{x', t) , x' = S{t) - \x\ 

where C(r) is given by eq. ( |6^ ) and initially d{x',0) « 
x'6{x'). Assuming that 

\dd/dx'\ > |d/T| and \ddldx'\ > \ddldT\ (103) 



and remembering eq. (|61| ) as the condition of cancellation 
of the dd/dx' terms, we find: 



Teg- 



d^d 



(104) 



dx'dr '■^ dx'^ 
Integrating over dx' from — oo to a given x' , we arrive at 
the classical diffusion equation 

dd 

The soluti on of this equation satisfies the second condi- 
tion in eq. (103) provided Tc is much less than the charac- 
teristic scale over which d{x' ,t) varies in x' . I ndiv idual 
Fourier modes of any function satisfying eq. ( 104 ) are 
damped as 



d{k,T) 



Mr) 



d{k,0) 



where 



4M= / 9ir')r,ir')dr' , 



(105) 



(106) 



in agreement with eqs. 

By eq. (99b), small scale CDM dynamics is not affected 
by Silk damping. 
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VII. CONCLUSIONS 

We have shown how to reduce the Unearized cosmo- 
logical dynamics of gravitationally interacting species in 
the conformal Newtonian gauge to a system of coupled 
wave equations. These equations indicate that the hnear 
evolution of not only the density fluctuations but also of 
the corresponding gravitational potentials (metric per- 
turbations) proceeds causally and locally, as opposed to 
the instantaneous action at a distance seemingly implied 
by the Poisson equation. In the fluid approximation, 
a disturbance in the gravitational potential propagates 
through the photon-baryon fluid at the speed of sound. 
The locality of the linear dynamics permits the efficient 
analysis of the perturbation evolution using the Green's 
function method. The Green's functions are simply the 
Fourier transforms of the familiar transfer functions, but 
the causal nature of the Green's functions provides in- 
sights that were previously unknown. 

We find that the Green's functions of primordial 
isentropic perturbations prior to photon decoupling are 
sharp-edged acoustic waves expanding with the speed of 
sound in the photon-baryon plasma. As shown in Fig. 3, 
much of the integral weight of the photon density pertur- 
bation is localized at the acoustic wavefront of the cor- 
responding Green's function. When the finite mean free 
path to Thomson scattering is accounted for, the photon 
density perturbation is broadened to the width of the Silk 
damping length. The photon dynamics in the presence of 
Thomson scattering is described by the Boltzmann equa- 
tion, which we have analyzed to first order in the photon 
mean free path. The result is a diffusive damping cor- 
rections to the radiation wave equation. We show how 
photon polarization affects Silk damping and provide an 
alternative derivation of the results of Kaiser, Ref. 

Another insight obtained with the Green's function 
method concerns the effects of baryons. The baryonic 
component of the plasma is responsible for a distinctive 
central dip in the Green's function for the gravitation- 



ally redshifted CMB photon temperature perturbation 
j6r + (p as shown in Fig. 5. This feature gives rise, 
Ref. , to the known effect of odd peak enhancement and 
even peak suppression, Ref. 0, in the CMB power spec- 
trum. It is readily understood in our analysis through 
the gravitational effect of the cold dark matter acting on 
the baryons. 

In this paper we considered perturbations only in pho- 
tons, baryons, and CDM. The Green's function method 
can be applied to the linearized phase space dynamics 
of an arbitrary number of species with rather general in- 
ternal dynamics and mutual interactions, including neu- 
trinos, Ref. jsj]. The advantage of the position space 
approach over the traditional Fourier space expansion is 
its simple and explicit treatment of acoustic and trans- 
fer phenomena underlying the dynamics of all particle 
species. In addition to providing an intuitive, compact 
framework for the many effects shaping the fluctuations 
of matter and radiation, this approach can give simpler 
analytical and faster computational methods for the still 
challenging aspects of cosmological perturbations, such 
as those of relic neutrinos. Application of the simi- 
lar wave equations and the Green's function approach 
to other areas of astrophysics where gravity influences 
acoustic or transfer phenomena might be fruitful. 



APPENDIX A: BOLTZMANN HIERARCHY 
IN POSITION SPACE 

Two independent physical potentials are needed to 
specify scalar metric perturbations in the general case, 
see e.g. any of Refs. g |, |, |, §. In the conformal 
Newtonian gauge, these are 0(r, r) and ijj{r, r) defined 
by eq. (^. The linearized Einstein equations in a flat 
universe reduce to the following equations, following cor- 
respondingly from the 0-0, 0-i, summed i-i, and trace- 
less i-j components of G^'^ = SnCTt"^, Ref. [|l|: 



2 2 _ (£) 



= 'i7rGa'^J2aPaSa , 

= 'iTrGa'^Y.aiPa + Pa)Ua , 



(Al) 



The energy density enhancement 6a and the velocity po- 
tential Ua of each of the matter or radiation species "a" 
are defined, as previously, in terms of the corresponding 
energy-momentum tensor components by eq. ([To|). The 
variables Spa and tTq , giving the isotropic and anisotropic 



components of stress perturbation, are defined by 

- S){Pa+Spa) + 

+ (Pa+Pa)^(v%-i<5i-V2)7r, (A2) 

where V* = V^, assuming negligible 3-space curvature. 
The anisotropic stress potentials tt^ vanish for perfect flu- 
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ids and we see from the last of eqs. ( |Al[ ) that the gravi- 
tational potentials (j) and ip are then equal. 

Six variables specify the coordinates of a particle in 
phase space at a given time. For them, we take the co- 
moving coordinates of the particle and the comoving 
momenta 



(A3) 



where pi are the proper momenta measured by a co- 
moving observer, Refs. The momentum coordi- 
nates Qi are not canonically conjugate to the variables r'; 
the canonical momenta of a particle of mass rua are 
Pi = rUadxi/V-ds'^ = (1 - ^)Qi- 

The particle density in phase space is specified by the 
canonical phase space distribution /a(?'*, Pj, t): 



dNa^ fa{r\Pj,T)d^r'd^Pj 



(A4) 



for every species of particles and their states of polariza- 
tion a. Time evolution of the phase space distributions 
is given by Boltzmann equation 



dfa 

dr 



dfa ^. dfa 

dr^ dq 



dfa 



dfa 

dr 



(A5) 



where fa is considered as a function of the coordinates , 
q=\qi\, rii = qi/q, and r. The energy- momentum ten- 
sor of the species a is given in the conformal Newtonian 
gauge by the following simple expression up to the first 
order of cosmological perturbation theory, Refs. H, |3^ : 



T!f,.= I d'p^P^fa 



with p° 



'Po = \/((?/a)2 + ml and p' 



Pi 



(A6) 



Qi/a. We 



later omit the species label a if referring to any sort of 
particles in general. 

As in the tight-coupling case, an arbitrary perturba- 
tion (5/(r, q, n, r) of the phase space distribution about 
the unperturbed distribution /o (q, r) can be expanded 
over plane waves, in which the value of Sf at a given q, 
n, and r remains constant along two spatial directions. 
To describe the dynamics of such a plane wave, we set the 
spatial coordinates y and z to vary along these invariant 
directions, and the coordinate x to vary along the remain- 
ing direction. The perturbation 5f(x,q,h,T) can be ex- 
panded over Green's functions satisfying the delta func- 
tion initial conditions ^/^^^ (x, g, n, r) — > A{q,n) 5d{x — 
xq) as t ^ with various initial locations xq and func- 
tional forms of A{q,n). In the translationally invariant 
background, the time evolution of the Green's functions 
is independent of xq. 

If ip is the angle in the y-z plane between the vec- 



tor {n. 



■y, ' 



and the axis y, we may expand an arbitrary 



A{q, n) in the Fourier series 



A{q,n^,ny,n^) = 



+ OC 



(A7) 



By definition, the scalar perturbations have initially zero 
A„i for m ^ 0. For a homogeneous and isotropic back- 
ground, subsequent evolution according to the Boltz- 
mann equation ( |A5| ) in the linear regime can not ex- 
cite m ^ components of 5f through scalar pertur- 
bations alone. Therefore, we will consider only the 
perturbations that are axially symmetric about x axis: 
5f = Sf{x,q,n^,T). 

The number of coordinates needed to describe Sf can 
be reduced further when the particle mass is zero and the 
relevant scattering cross-sections are energy independent. 
Then one can work with the energy averaged perturba- 
tion of the distribution, Ref. ISgI , 



r., X J q'^dqq5f{x,q,n^,T 



j q^dqqfoiq) 

OO 

^Y.^2l + l)Fi{x,T)Pi{n.,) , 

1=0 



(A8) 



where Pi is a Legendre polynomial. The energy density 
enhancement, mean velocity, and anisotropic pressure of 
the particles described by an energy averaged distribu- 
tion F are given by 



S 

-Vit 

V^TT 



j\^Fin,)^Fo , 
lj\^n^F{n,)^^^F, 



(A9) 



for every massless species a. In these and the follow- 
ing equations V = since we are considering plane- 
parallel perturbations that are constant in the y and 
z dire ction s. Eqs. (A9) follow direct ly f rom the defini- 
tions (|lOa|, P^,|A^) and the formula ( |A6| ). 

Continuing our practice of reducing the number of gra- 
dients in the evolution equations, we define for all I 



F = {^l)'Vfi 



(AlO) 



The evolution equations for the variables in eqs. (|A9|- 
AlO) can be d erived in position space from Boltzmann 
equation (A5), Ref. or written down by replacing 
momenta k by spatial gradients in previously derived 
Fourier space equations, e.g. Ref. |6|. We give the cor- 
responding hierarchy equations in the conformal Newto- 
nian gauge (H) for neutrinos, photons, baryons, and CDM 
particles. 

The evolution of neutrinos is given by 



■ 5^ + V^TT^ 



(All) 



fvl 



I 



21 + 1 



fu(i- 



(i-i) 



l + l 
21 + 1 

3 f 

4 J fl 



I > 2 



with TT^ = -i f,^2 and v 

The collision terms for Thomson scattering of photons 
and baryons were derived in Refs. [p7|. The resulting 
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equations read 

= ^ V^u^ + 4?/' , 

= - S-f + V^TT^ + </> (u-y — Ub) , (A-12) 

4 Tc 

Again, TT-y — \ f^2, = I/7I7 and Tc is defined by 
eq. ([75I). The quantities gji{x,T) and g^i{x,T) are de- 
fined as 



G70 — V^5^o ) 

G^i = — V^5^i , 

G^i = (-l)'V'5^z , l>2 



(A13) 



where G is the energy averaged distribution describing 
the difference of two linear photon polarization compo- 
nents. The equations of their evolution are 



1 



1 



(A14) 



57O = V - — |g^o - 2 + .970 + 372) } , 
^ _ I ^ 2 1 , 

^71 — 2 ff70 + 2 572 - — 57I 

I ^ + 1 w2 

^ 27TT'^"('-^) + 27TT^ ^''^"-'^ " 

~ ~" (ff^i - ^ (/72 +57O +572)1 , ^>2 



where g^i = V^^^i in the last equation at / = 2. 

Applying momentum conservation to photon-baryon 
interactions, one can easily modify the non-relativistic 
fluid equations of baryon evolution to include Thomson 
scattering: 



(A15) 



Ub = Ub 

a 



1 4 

Tc 3/3y 



{ub - Uy) 



The parameters y and /? are defined by eq. (|T^). Aside 
from the distinction between and ■0 metric perturba- 
tions, the equations of CDM evolution are the same as 
in the fluid model: 



Uc = Uc 

a 



-30 



(A16) 
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